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We consider the weakly first order phase transition between the isotropic and ordered phases of 
nematics in terms of the behavior of topological line defects. Analytical and Monte Carlo results are 
presented for a new coarse-grained lattice theory of nematics which incorporates nematic inversion 
symmetry as a local gauge invariance. The nematic-isotropic transition becomes more weakly first 
order as disclination core energy is increased, eventually splitting into two continuous transitions 
involving the unbinding and condensation of defects, respectively. These transitions are shown to 
be in the Ising and Heisenberg universality classes. A novel isotropic phase with topological order 
occurs between them. 



I. INTRODUCTION 

The modern theory of critical phenomena [1] empha- 
sizes the role of symmetry and dimensionality in de- 
termining the long-wavelength behavior at phase tran- 
sitions. From this point of view, the properties of the 
ferromagnet-to-paramagnet [1] and nematic-to-isotropic 
[2] order-disorder phase transitions might be expected to 
be quite similar. In both cases, a continuous rotational 
symmetry is broken, leading to the same sort of long- 
wavelength Goldstone modes, spin waves and director 
modes, respectively. Because of this similarity, an ex- 
pansion [3] in e = d — 2 for spatial dimensions d near two 
predicts that both transitions are continuous and have 
identical exponents to all orders in e. 

Prior to our earlier work [4], it was generally believed 
that this prediction was wrong in three dimensions, where 
the ferromagnet to paramagnet transition is continuous, 
whereas all observed nematic to isotropic transitions are 
first order, albeit usually only weakly so. 

These empirical facts are in accord with the predic- 
tions of Landau theory [2], although that treatment fails 
to explain the ubiquitous weakness of the first order ne- 
matic transition. The Landau approach treats these two 
order-disorder transitions using qualitatively different or- 
der parameters, and apparently disregards the inherent 
similarities of the ordered states. In this paper we will 
focus on the role of topological defects in the nematic 
phase transition, since "disclinations" - the line defects 
which occur in nematics but are not allowed in magnetic 
systems - are the sole topological distinction between the 
two kinds of ordering. 

This paper presents a thorough analysis of a new model 
[4] of the three-dimensional nematic-isotropic (N/1) tran- 
sition in which a pivotal role is played by disclinations. 



Our model allows us to independently vary the local ne- 
matic stiffness and the disclination core energy. When 
the defect core energy is large, our nematic contains only 
a few, thermally activated, small disclination loops, and 
is therefore similar to the magnet, where such loops are 
topologically forbidden. 

We find that sufficiently strong suppression of de- 
fects leads to behavior that contradicts the predictions 
of Landau theory. In particular, the first order nematic- 
isotropic transition splits into two continuous transitions, 
with a novel intervening "topologically ordered" isotropic 
phase (see figure 1). The transition from the conventional 
isotropic phase (denoted by '/') into the new topolog- 
ically ordered isotropic phase (T') belongs to the uni- 
versality class of the three-dimensional Ising model, and 
that from the topologically ordered phase to the nemat- 
ically ordered phase ('N' for "nematic") to the three- 
dimensional, three-component (n = 3) Heisenberg uni- 
versality class. Recent experiments by Roux et al. [5] 
may have revealed such a phase diagram. 

The organization of this paper is as follows. In sec- 
tion II we further discuss the the conflicting approaches 
to the N/I transition based on Landau theory and 2 + e 
expansion, and describe the disclinations which are the 
sole topological distinction between ferromagnets and ne- 
matics. Motivated by these considerations, we present in 
section III two models which explicitly suppress these de- 
fects, with the purpose of making the nematic more like 
the ferromagnet. The phase diagram of this second model 
(figure 3), or more specifically the existence of three dis- 
tinct phases and the natures of the transitions between 
them, is determined in sections IV- VII. Analysis of three 
limits which reduce to already understood models is pre- 
sented in section IV, followed, in section V by a character- 
ization of the three phases in terms of defect line tension, 
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and the introduction of order parameters in section VI. 
The robustness of the exact limits analyzed in section IV 
is demonstrated in section VII, thereby establishing in 
particular that the nematic to physically isotropic (but 
topologically ordered) transition remains continuous and 
in the Heisenberg universality class over a non-vanishing 
range of parameters. This section is fairly dense and can 
be skipped on first reading. Section VIII shows that per- 
turbations such as space-spin coupling that are absent 
from our model yet present in real nematics do not alter 
our conclusions. Our Monte Carlo results on the new ne- 
matic gauge theory, including a finite size scaling analysis 
of the phase transitions, are presented in section IX, and 
corroborate the analysis of earlier sections. A separate 
paper [6], Paper II, presents detailed predictions for crit- 
ical phenomena at the new transitions predicted here. 



II. SYMMETRY, EPSILON EXPANSION, AND 
LANDAU THEORY 

Fluctuations. Ferromagnets and nematics both 
spontaneously break a continuous rotation symmetry: 
the elementary moments of a ferromagnet are preferen- 
tially aligned parallel to a common vector in "spin space," 
while in a nematic the long axes of the constituent mole- 
cules preferentially align with a common axis in physical 
space. By "axis" we mean a headless vector, usually 
called a "director." Note that there are two directions 
associated with an axis, so the nematic has a local "up- 
down" symmetry which is lacking in the ferromagnet. 

The "order parameter space" of an ordered medium 
contains the set of "directions" that are associated with 
the broken symmetry. The order parameter space for 
the Heisenberg magnet is the unit sphere S 2 . For the ne- 
matic, the order parameter space is the projective plane 
IP 2 , which is simply the unit sphere with antipodal 
points identified. There is evidently a two-to-one map- 
ping from S 2 to BP 2 which is a local isometry. Thus small 
fluctuations of both order parameters therefore have sim- 
ilar phase space measures. 

In an ordered phase, the effective free energy for small 
fluctuations can be written in terms of gradients of the 
order parameter. Since the nematic and ferromagnetic 
order parameter spaces are locally isometric, the asso- 
ciated long-wavelength physical properties are also the 
same. Indeed, the Frank free energy density for nemat- 
ics, 

■^"gradW = i[#i(V • n) 2 + K 2 (n ■ (V x n)) 2 

+ K 3 (n x (V x n)) 2 ], (1) 

is an acceptable spin- wave Hamiltonian for a ferromagnet 
with space-spin coupling. (If all of the Ki are equal, 
eq. (1) becomes simply (A'/2)(3 ,n ( g)(3"n /3 ), which is 
invariant under independent rotations of space and of 
the order parameter.) 




FIG. 1. The projective plane MP 2 . Points on the equa- 
tor and in the 'missing' hemisphere are identified with their 
antipodes, so that it does not have a boundary, contrary to ca- 
sual appearances. A closed, but non-contractible, path start- 
ing and ending at point x is depicted. 

An immediate corollary of this discussion is that any 
treatment of nematic ordering which considers only non- 
singular fluctuations about a uniform, nematically or- 
dered state will necessarily be identical to the correspond- 
ing treatment of a ferromagnet. This is precisely the ap- 
proach of the e expansion [7] . The existence of a critical 
fixed point of order e implies that both transitions can 
[8] be continuous, and would share the same universality 
class. 

Topology. Although the local fluctuations of the or- 
dered nematic and ferromagnetic states are quite similar, 
certain configurations of the two systems are very differ- 
ent. Suppose for example that we arbitrarily assign an 
arrowhead to the director at some point in a nematic, and 
try to extend this smoothly to a continuous vector field 
that is consistent with the given director field. For small 
fluctuations this is clearly possible; traversing a closed 
path in physical space maps out a closed path on the 
unit sphere S 2 . Such configurations are "homotopically 
trivial." Note that all smooth magnetic configurations 
are trivial in this sense. Homotopically trivial nematic 
configurations can therefore be placed in correspondence 
with related magnetic configurations. Since their ener- 
gies will be both be governed by the elastic energies of 
eq. (1), we expect these configurations to make similar 
contributions to the partition function of their respective 
systems. 

The nematic configurations which cannot be related 
to magnetic configurations in this manner are those for 
which an ambiguity in sign arises at the completion of a 
closed path in physical space. That is, the image in MP 2 
of a closed path in the physical system will be closed if 
the director field is continuous, but its "lift" to S 2 travels 
from one point to its antipode and is not closed. A closed 
curve in the nematic that corresponds to a homotopically 
nontrivial loop in HP 2 necessarily encircles a singularity 
of the order parameter, since shrinking the path in real 
space induces a deformation of the corresponding path in 
the order parameter space. Such a configuration is said 
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(b) (c) 
2. (a) A cross section through a disclination. When 



the defect is encircled, the local molecular axis rotates through 
180°. On a lattice the defect can occur inside a plaquette. The 
presence (b) or absence (c) of a defect depends on the product 
of the Uij around the plaquette. 



to contain a topological line defect known as a "disclina- 
tion." The topological character of these defects endows 
them with considerable dynamical stability. This essen- 
tial distinction between nematic and magnet is quantified 
by the fundamental group The fundamental group 
7Ti(5 2 ) of the sphere is the trivial group with only an 
identity element, whereas that of the projective plane, 
7Ti(ffiP 2 ), is the group Z 2 (integers mod 2 under addi- 
tion). 

A cross-section of a nematic configuration with a discli- 
nation is depicted in figure 2. Related configurations in 
the magnet have very large energies, because there are 
180 degree discontinuities in the corresponding vector 
field. Thus configurations containing defects will make 
distinct contributions to the partition functions of the 
nematic and magnet, and provide a possible explanation 
for the fundamental difference between the two order- 
disorder transitions. These disclinations are very impor- 
tant for the nature of the nematic to isotropic transi- 
tion. That this is so is suggested by earlier work (for a 
very readable review see reference [9]; also see references 
[10-12]) showing that many phase transitions proceed by 
an unbinding of topological defects. 

The most famous example is the explanation of the 
transition in the two-dimensional XY model as an un- 
binding of point vortices by Kosterlitz and Thouless. 
[13,14] In three dimensions, both the XY model and the 
superconductor can be cast as problems of interacting 
vortex loops [15]. This produced the discovery that the 
superconducting-normal transition can be continuous (as 
experimentalists are well aware!), in contrast to the pre- 
diction of the 4 — e expansion that fluctuations of the 



electromagnetic field drive a first order [16] transition. 
Monte Carlo studies [17] have found that suppression of 
these loops (by adding a sufficiently large core energy) 
constrains the model to remain in its ordered phase, and 
prevents a phase transition to a disordered state. 

In the three-dimensional Heisenberg model, there are 
no stable line defects, but point defects ("hedgehogs") are 
possible. Suppression of these point defects in the three- 
dimensional Heisenberg model has been studied in simu- 
lations by Lau and Dasgupta [18], in which the transition 
was destroyed by a sufficiently large energetic penalty for 
hedgehogs. Fucito and Solomon [19] likewise found that 
the three dimensional XY transition was destroyed by 
suppressing vortex lines. 

Landau theory. In Landau theory, the free energy is 
expanded in powers of an "order parameter [21]." Since 
the nematic has global inversion symmetry its order pa- 
rameter cannot be a vector: there is no macroscopic di- 
rection selected by the ordered state. There is, however, 
a preferred axis, which leads to anisotropy in e.g. mag- 
netic and dielectric susceptibility. The order parameter 
can then be chosen to be a traceless symmetric tensor 
Q a f3- (For this discussion of Landau theory we restrict 
ourselves to three-dimensional systems, so a and j3 here 
run from 1 to 3.) . For weak nematic order Q a p is propor- 
tional to the anisotropy of the susceptibilities and light 
scattering, and can be thought of as the quadrupole mo- 
ment of the local distribution of molecular orientations. 

The Landau theory of the nematic-isotropic transition 
is readily constructed in terms of Q by writing down all 
rotationally invariant scalars. To fourth order, the Lan- 
dau free energy is 

^"lg = a(T - T ) Tr [Q 2 ] - bTr [Q 3 ] + 2c Tr [Q 4 ]. (2) 

(Note that the fourth order term is unique, since 
(Tr Q 2 ) 2 = 2 Tr Q 4 for any traceless symmetric three- by- 
three matrix.) Q has five independent components, and 
the quadratic term is simply the sum of their squares. 
Thus, if b = 0, a system governed by this free energy 
will undergo a second-order transition at To in the 0(5) 
universality class. The cubic term is relevant at the asso- 
ciated fixed point, however, and there certainly appears 
to be no reason for b to vanish generically. When 6^0 
the system undergoes a first order transition at 
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For a magnet, however, the corresponding free energy 
cannot have a cubic term, since inverting the magneti- 
zation is a symmetry of the system and therefore cannot 
change the free energy. 

For b 0, the free energy (2) is minimized by a uniaxial 
tensor 



Qa/3 = S{n a np - l/36 a p), 



(4) 



where n is a unit director. The free energy then reduces 
to 
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Uniaxiality is stable against the introduction of higher 
order terms in eqn. (2) for sufficiently weak nematic or- 
der. 

Nematic transitions observed in nature are first order 
(score one for Landau), but generally only weakly so [2]. 
This is revealed by light scattering, which shows nearly- 
critical fluctuations, and means that the transition tem- 
perature T* is close to the spinodal To. The jump in the 
order parameter at the transition, 
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is not necessarily small. The Landau theory gives no 
clue as to why the transition should generally be so weak 
despite significant variation in microscopic properties, 
though it has been suggested [2] that the aspect ratio of 
the constituent molecules may provide a natural "small 
parameter." 

Motivation for our work. We have seen that the 
e expansion, which ignores disclinations, predicts a con- 
tinuous ncmatic-to-isotropic transition, while the Lan- 
dau theory, which implicitly includes disclinations (since 
it averages, albeit crudely, over all configurations of the 
full nematic order parameter, including those with de- 
fects), predicts a first-order transition. Is it, then, the 
disclinations that drive the nematic-to-isotropic transi- 
tion first-order? And if so, can the continuity of the 
transition be restored by suppressing these defects? Our 
answer to both questions is a resounding "yes!" These 
conclusions are reached by considering a new theory of 
nematics which incorporates both spin-wave and topo- 
logical fluctuations on an equal footing. We find that 
the strength of the first-order nematic-to-isotropic tran- 
sition is reduced as the defect core energy is increased. 
For large enough defect suppression the transition splits 
into a pair of continuous transitions, with a qualitatively 
new intervening phase (see figure 3). This new phase is 
isotropic ({Qap) = 0), like the fully disordered phase, but 
possesses distinct topological characteristics. 

Experimental verification. To test these ideas ex- 
perimentally one must have independent control over the 
disclination core energy, corresponding to K, and the lo- 
cal nematic interaction strength corresponding to J. This 
is likely to be difficult, since in real materials both mi- 
croscopic energies arise from the same interactions and 
entropic configurations. Nevertheless one can imagine a 
"defectophilic" impurity which would accumulate in de- 
fect cores, and thus increase their entropic cost. Sim- 
ilarly, long "defectophobic" molecules might align with 
the nematogens and make formation of defects more dif- 
ficult. Observation of weakening (strengthening) of the 
first order transition as defects are suppressed (favored) 
would provide partial support for our scenario. 

Another approach is to add objects to the nematic 
which disorder it without favoring the creation of discli- 
nations. In this way, the nematic order is destroyed at 
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FIG. 3. Phase diagram for the Z2 gauge 
plus three-component spin lattice model. The phase bound- 
aries are very straight except near their intersection. The two 
nematically disordered phases are not distinguishable by any 
local quantity. The dotted line denotes a first-order transition 
and the unbroken lines continuous transitions. 



a lower temperature, where disclinations are effectively 
suppressed. (That is, Ec/UbT increases while the defect 
core energy E c remains the same.) This effect may ex- 
plain the experiments of Roux et al. [5], who added small 
polystyrene spheres to a nematic. These spheres favor the 
creation of point defects (hedgehogs), and therefore tend 
to destroy nematic order. Roux et al. appear to observe 
critical opalescence, suggesting a continuous transition 
between two isotropic phases for this system; it is tempt- 
ing to identify this transition with the 1/ T transition in 
figure 3. Detailed experimental tests of the critical prop- 
erties described in Paper II would confirm this identifi- 
cation. 

Unfortunately, the new topologically ordered phase we 
predict will be difficult to observe directly because it 
is isotropic, and therefore quite similar to the conven- 
tional isotropic state. Nevertheless, the passage into and 
out of this new phase can be recognized by the associ- 
ated critical singularities. These singularities are deter- 
mined by the universality classes of the transitions, which 
are three-dimensional Ising for the transition between 
the fully disordered and topologically ordered phases 
(7/ T) and three-dimensional Heisenberg for the transi- 
tion between topologically and nematically ordered states 
(T/N). Detailed predictions can be found in Paper II. 



III. LATTICE MODELS FOR NEMATIC 
SYSTEMS 

In this section we introduce two lattice models for ne- 
matic media which explicitly include energetic penalties 
for disclinations: a "spin-only" model and a lattice gauge 
theory. The behavior of the spin-only model is less in- 
teresting, but it clarifies the motivation of for the lattice 
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gauge theory. Monte-Carlo studies of this latter model 
is reported in section IX. Analytic results for the lattice 
gauge theory will be studied in more detail in sections IV 
and VII. 



A. Modified Lebwohl-Lasher Model 

The simplest lattice model of the nematic-isotropic 
transition was introduced by Lebwohl and Lasher [22], 
and is based on the continuum Maier-Saupe model [2]. 
The local nematic degrees of freedom are represented by 
unit vectors Si that are assigned to the sites of a regu- 
lar lattice (cubic for convenience). The free energy of a 
configuration of such spins is given by 

f3H N = -Jjv^Si-S,-) 2 (7) 

where (i, j) denotes nearest-neighbor pairs of lattice sites. 
Since the spins occur quadratically, there is ZocaHnversion 
symmetry, as appropriate for a nematic model. That is, 
the spin at site i may be negated without changing the 
energy of the system. This local symmetry foreshadows 
the introduction of a gauge description of nematics (see 
below). 

To consider topological defects, we introduce a defect 
counting operator 

Vijki = ijl - sgn[(S, • S j )(S j • S k )(S k • S ( )(S ( • SO] J, 

(8) 

which is unity if a defect pierces the plaquette (ijkl), and 
zero otherwise. This four-spin operator evidently retains 
the local inversion symmetry of eqn. (7). 

If the term in square brackets vanishes, then by a 
process of flipping spins (to which the nematic config- 
uration, hence the energy, is insensitive), we can make 
all the factors (Si • Sj) positive. Then the smooth in- 
terpolation which takes the shortest route between the 
corresponding values for the director (on MP 2 ) is homo- 
topically trivial. Thus, there is no disclination threading 
the plaquette (ijkl). In contrast, in the presence of a de- 
fect there is always one leftover negative factor, and the 
defect number is unity. The associated path on MP 2 is 
nontrivial. 

A defect core energy is then 

(3H D ee K D J2^ijki (9) 
□ 

where the sum extends over all elementary plaquettes 
□ = (ijkl). A Monte Carlo study shows that the tran- 
sition weakens (as measured by the order parameter dis- 
continuity) as Kd increases. There are only two states 
(nematic and isotropic) in the phase diagram of the 
model of eqns. (7-9), and the transition between them 
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FIG. 4. The link variables Uij provide information on the 
rotation of the local molecular axes between sites i and j. 

is always first order, as expected by Landau theory. For 
large enough Kjj, however, the system remains ordered 
even for vanishing bare nematic stiffness Jn- [23] Thus 
the elimination of defects leads to a long-wavelength 
renormalized nematic stiffness which is larger than the 
value required for long-range order to appear. This re- 
sult mimics related work on the three-dimensional XY 
[17] and Heisenberg [18] models, where elimination of de- 
fects leads to long-range order even in the absence of a 
bare order-parameter stiffness. This model, therefore, 
nowhere exhibits a continuous nematic-to-isotropic tran- 
sition. In the next subsection, we consider a different 
model, but in the same spirit as this one, which does 
achieve this goal. 

B. Lattice Gauge Theory 

The model of equations (7-9) lives on a lattice. To 
deal with topological defects, however, we had to con- 
sider what happens between lattice sites. To take this 
into account from the beginning, we envision the lattice 
as embedded in a (continuum) nematic fluid. A coarse- 
grained director at each lattice point can be constructed 
from the mean molecular axis within a region whose ra- 
dius is somewhat less than the lattice spacing. This pro- 
cedure can be carried out as long as the correlation length 
is larger than the lattice spacing. In our model this local 
molecular axis is associated with a vector Si at each site 
i, as in the Lebwohl-Lasher model. Once again, the sign 
of Sj at each site is arbitrary. 

Now consider the variation of the local molecular axes 
between two lattice sites i and j (see figure 4). Beginning 
at Si on the unit sphere and following the variation of the 
director, we trace out a curve which ends at either Sj or 
— Sj. Thus we arrive at two homotopy classes of paths 
as discussed in the introduction. Those in each class 
are deformable into one another (while keeping endpoints 
fixed). In our lattice description of the nematic fluid, we 
must retain this information regarding the director field 
between sites; it is needed to define topological defects, 
which require a notion of continuity. The presence or ab- 



5 



sence of a defect within a plaquette is determined by the 
homotopy class of the director field around it, which in 
turn is determined by splicing together the classes corre- 
sponding to the links making up the plaquette. 

The local degree of freedom associated with each link 
(ij) is represented by Uij, which can be either +1 or — 1. 
A value +1 of Uij means that one can assign a continu- 
ously varying orientation to the molecules between i and 
j which matches up with both Sj and Sj (fig. 4a). If the 
orientation can be chosen to match either S t or Sj, but 
not both, Uij is —1 (fig. 4b). 

This approach differs crucially from the modified 
Lebwohl-Lasher model of eqns. (7-8), which implicitly 
assumes that the local nematic axis follows the shortest 
route between its values on the lattice sites. If the lattice 
is not too coarse, this is justifiable by energetic consider- 
ations which suggest that these configurations will dom- 
inate the partition function. When defects are strongly 
suppressed, however, this need no longer be the case. 
The extra energy cost incurred by "taking the long way" 
may be more than compensated by avoiding the imposed 
defect core energy. Thus in a strongly defect-suppressed 
regime, we must keep track of configurations in both the 
defective and defect-free classes. 

By introducing the link variables Uij, the local Z 2 in- 
variance of a nematic fluid is expressed as a local gauge 
invariance. Choosing fa = ±1 independently for each 
site i, the energy of any configuration is unaltered if the 
simultaneous transformations 

Si ► </>iSj 

Uj -► OiOjUij (10) 

are performed on the site and link variables. The simplest 
Hamiltonian including a defect-suppression term which 
respects this local gauge symmetry is 

(in = -jJ^UijSi ■ Sj - KY.UijUjkUuUH, (11) 

m n 

where the second sum is over all elementary plaquettes 
□ = (ijkl). The partition function is found from e _/3W 
by integrating over spins {Si} and summing over link 
variables {Uij}. 

The first term in eqn. (11) is a nematic interaction 
which favors minimal variation of the director along link 
(ij). For example, the configurations depicted in figs. 3a 
and 3c have lower energy than those in 3b and 3d. 

The second term is a defect core energy analogous to 
cqn. (8). If the product of the link variables Uj around 
a plaquette is +1, then there is a smooth pattern of ori- 
ented molecular axes along the links which does not use 
the head-to-tail symmetry of the nematic. This pattern 
can be continuously extended to the entire interior of the 
plaquette, indicating the absence of a disclination. If the 
product of the link variables Uj around a plaquette is 
— 1, however, the local molecular axis rotates by 180° as 
the plaquette is encircled (compare figures 2b and 2c). 



This is most easily seen by using a gauge transformation 
to set as many of the links on the plaquette to +1 as 
possible (all of them if the product of links is +1, all but 
one otherwise). 

How do the defect suppression terms in equations (8) 
and (11) compare? From the first term in eqn. (11), we 
see that it is energetically preferable for sgn(S; • Sj) to 
be equal to Uij. Thus UijUjkUkiUu (sites i,j,k,l around 
a plaquette) is a rough measure of Y[ sgn[(S; • Sj)]. This 
equivalence is stronger for larger J, since then S is less 
likely to take the long way between its values on the lat- 
tice sites. Later we will see that the resemblance between 
the gauge model and that of equations (7,8) is greatest 
in the limit of small K. We have performed Monte Carlo 
simulations of the two related lattice nematic models of 
eqns. (7-9) and (11), and our results are discussed in 
detail in section IX. 



IV. EXACT EQUIVALENCES 

In this section, we show that the gauge model must 
have three distinct phase transitions at the points M, I 
and H on the phase diagram of figure 3. This provides 
strong preliminary evidence that there are at least the 
three phases which actually appear there. 



A. Complete defect suppression: K = oo 

When K = oo, only defect-free configurations con- 
tribute to the partition function. The product of the 
Uij around each elementary plaquette, hence around any 
closed path, then must be +1. The defect density is zero 
everywhere, and the gauge can be chosen [24] so that 
Uij = +1 for each link (ij). Specifically, this can be done 
as follows: For each site i, pick a path of links P(i) from 
the origin to i and define 

o-i= n Uki - ( i2 ) 

(ki)eP(t) 

Since the product of link variables around any closed path 
is +1 at infinite K , <7; so defined is independent of the 
specific choice of path P(i). Furthermore, 

Uij = cr t o-j, (13) 

as can be seen by constructing P(i) as the concatenation 
of P(j) with the single link (ij). Thus, the partition 
function becomes 

Z=J2 [{dS}ex P {-jJ2(^)-(°jSj)} , (14) 

and the tr's embody the freedom of gauge choice. Chang- 
ing to new variables Sj = er^Si, each term is independent 
of the ct's and the partition function is 2^ times that 



6 



for an n = d = 3 Heisenberg model. Thus for K = oo 
there is a second order transition in the three-dimensional 
Heisenberg universality class at [20] J ~ 0.693. 

We show in the next section that this Heisenberg tran- 
sition persists for large but finite K. Finite size scaling 
analysis of our Monte-Carlo simulations confirms that the 
transition remains in the Heisenberg universality class 
out to the multicritical point M. Thus "defect fugacity" 
is irrelevant, in the renormalization group sense. 

For complete defect suppression (K = oo), in the gauge 
with all Uij equal to +1, the ordered phase at large J is 
characterized by a non-zero total magnetization. One 
may instead check for a non-zero limit of the usual spin- 
spin correlation function lim|j_j|^ 00 (Si • Sj) > 0. This is 
generalized to arbitrary gauge in the next section. 

B. Pure gauge theory: ,7 = 

For J = the spins decouple completely, and our 
model (eqn. 11) reduces to the three-dimensional pure 
Ising lattice gauge theory on a cubic lattice. This 
gauge theory can be mapped onto the ordinary three- 
dimensional Ising spin model by a duality transforma- 
tion [25]. The critical coupling J c = 0.22165 ± 0.00003 
for the three dimensional Ising spin model (from Monte 
Carlo [26] and high-temperature series [27] studies) then 
implies a critical point in the Ising universality class at 
K ~ 0.765. While J = 0, K ^ is of course an unphysi- 
cal limit of our model, the value of the mapping is in its 
use as a starting point from which we will argue that the 
transition persists and remains in the Ising universality 
class, for small J > 0. Finite size scaling studies of our 
Monte-Carlo simulations indicate that this Ising transi- 
tion persists all the way to the multicritical point M of 
figure 3. 

C. No defect suppression: K = 

When K = 0, the Ising link variables on different links 
are decoupled from one another. It is then trivial to 
trace over Uj = ±1 and obtain an effective Hamiltonian 
for the spins alone (subtracting a constant TV log 2): 

#eff [Si] = ln { Y, e MJUij Si-Sj)\-N log 2 

= 5^1n cosh(JS i -S j ) 

(ij) 

= 5]ij 2 (S l -S J ) 2 + 0(J 4 ). (15) 

(ij) 

Local gauge invariance of eqn. (11) guarantees that trac- 
ing over the £/'s will generate an effective Hamiltonian for 
the spins which is even in each S separately, i.e. which 
has local inversion symmetry. 



Like the Hamiltonian (eqn. 7) of the Lebwohl-Lasher 
model, f/eff is a function of (Si • Sj) 2 ; the resemblance is 
even stronger when it is Taylor expanded in J. Thus it is 
not surprising that, also like the Lebwohl-Lasher model, 
it has a single, first order phase transition between a 
nematic and an isotropic phase, indicated by the point N 
in figure 3 at J « 1.9. 



V. PHASES OF THE NEMATIC GAUGE 
THEORY 

In this section, we describe the macroscopic distinc- 
tion among the three phases depicted in figure 3. In the 
nematically ordered phase N, rotational symmetry break- 
ing and nonzero order parameter are accompanied by a 
non-zero helicity modulus [28], or spin- wave stiffness. In 
our model, the one-Frank-constant approximation is ex- 
act. In a more general model than ours, which included 
space-spin coupling, the single spin-wave stiffness is re- 
placed by the three Frank constants. This means that if 
the nematic fluid is confined to a box of side L, imposi- 
tion of boundary conditions with a relative angle 9 < ir 
between the directors on opposite faces of the box raises 
the free energy density by T(9/L) 2 over its value for pe- 
riodic boundary conditions (i.e., 9 = 0). The helicity 
modulus T is positive in the ordered phase and increases 
with the degree of ordering. Note that for fixed 9, the 
difference in free energy density vanishes in the thermo- 
dynamic limit L — ► oo. In the one-Frank-constant ap- 
proximation, T = K/2, where K is the Frank constant. 

An analogous measure of the free energy cost of discli- 
nation lines can also be constructed via response to 
changes of boundary conditions. Imagine a cylinder of 
height L and radius R filled with a nematic fluid. Con- 
sider two boundary conditions: 

(a) No disclinations are permitted to pierce the bound- 
ary, so that all defects must form loops contained 
entirely within the cylinder. 

(b) A single disclination is forced through the centers 
of the upper and lower faces of the cylinder, but 
no other defect lines are allowed to pierce the sur- 
face. Internal disclination loops coexist with the 
externally imposed defect. 

The three phases appearing in figure 3 are distin- 
guished by the dependence of the free energy difference 
8T between boundary conditions (a) and (b) on the ra- 
dius and height of the cylinder. Note that, as with the 
manipulation of boundary conditions used to measure 
the helicity modulus, this free energy difference is not 
extensive. 

Nematic phase. In the nematically ordered phase N, 
the extra defect imposes a variation of the director arbi- 
trarily far from the core, since the director must undergo 
a net rotation of ir along any path encircling the defect. 
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The free energy difference between the two boundary con- 
ditions is therefore governed by the helicity modulus, viz., 

= C N L\n(R/a), (16) 

where Cjv( J, K) = irT/2 is the long wavelength nematic 
stiffness and a is the defect core size. The defect line ten- 
sion SF/L is therefore logarithmically divergent with sys- 
tem size whenever long-range nematic order exists, and 
vice versa. 

A similar calculation shows that the interaction en- 
ergy between a pair of externally imposed defects has a 
logarithmic dependence on separation. Thermal fluctua- 
tions generate spontaneous defect loops within the cylin- 
der, but these loops have a strong energetic preference to 
avoid director variation far from their cores. They can 
achieve this by being small or binding in pairs. In this 
phase, the defects can be said to be confined. 

For small enough J (that is, J below a -/^-dependent 
threshold J C (K), the long- wavelength nematic stiffness 
vanishes and the free energy difference between (a) and 
(b) is no longer logarithmic with R. Depending on the de- 
gree of defect suppression, however, there are two distinct 
regimes which correspond to the two isotropic phases / 
and T of our model. 

Topologically ordered phase. Let us consider a 
state without long-range nematic order but with suf- 
ficient short-range order to permit the coarse-graining 
leading up to the gauge model of eqn. (11). If the defect 
core energy K is sufficiently large (e.g., K = oo), then 
the local molecular axes can be consistently assigned a 
continuous orientation throughout the system even in the 
absence of long-range order. This assignment effectively 
converts the short-range ordered nematic director field 
into a short-range-ordered vector field. We call this novel 
phase the topologically ordered state. Its physical prop- 
erties are isotropic (since it lacks long range order), but 
it is nevertheless distinct from the usual isotropic phase. 
This distinction is quantified by different asymptotic be- 
haviors of the defect free energy. 

Despite the absence of long-range nematic order in 
the large- K disordered phase T, the free energy cost per 
unit length of disclination remains nonzero. The bend- 
ing imposed on the spins by the presence of the defect is 
screened out over the correlation length, since the helic- 
ity modulus is zero. Near the defect core, however, the 
extremely tight bending (along with the bare defect core 
energy K) produces a non-vanishing defect line tension 
which is independent of the radius of the cylinder. The 
free energy difference between the two boundary condi- 
tions is given by 

6 Ft ~ G T L, (17) 

where Ct(J,K) is a constant which depends on the bare 
core energy and the elastic energy within a nematic cor- 
relation length of the core. 



In contrast with the nematic phase, in the topologically 
ordered state there is no long-range interaction between 
disclinations mediated by the nematic fluid, since long 
range nematic order is not present. Defects are there- 
fore unbound, and loops of finite extent will proliferate. 
Viewed from far away compared with their linear extent, 
these loops are either topologically trivial or the point 
defects of the nematic (hedgehogs). The picture of this 
phase as a gas of unbound hedgehogs is similar to that 
for the (d = n = 3) Heisenberg model in reference [18]. 

If the nematic interaction J is now increased, nematic 
order will develop from a state without topological de- 
fects. Since the relevant configurations are those of a vec- 
tor field, an ordering transition in the Heisenberg univer- 
sality class is expected. This is discussed more carefully 
in sections IV A and VII C. If the renormalized defect fu- 
gacity is small enough that arbitrarily large defect loops 
do not proliferate, the transition will be in the Heisen- 
berg class along the entire line MH, as confirmed below 
by Monte Carlo studies. 

Isotropic phase. As the core energy K is diminished 
at fixed, weak nematic interaction J, thermal fluctuations 
will cause the defect line imposed by boundary condition 
(b) to meander through the system. The system will 
then gain an entropy proportional to its length L, and 
the defect line tension will diminish. At a critical point 
the line tension will vanish, and the free energy differ- 
ence between (b) and (a) will become independent of the 
dimensions of the cylinder for large cylinders: 

T,=Cj(J,K). (18) 

We will show below that this transition can be under- 
stood as an Ising lattice gauge theory whose critical point 
lies in the universality class of the three-dimensional Ising 
model. 

In the resulting disordered phase /the nematic stiffness 
and the defect line tension both vanish, so the free energy 
difference between the two boundary conditions remains 
finite as L and R tend to infinity. In this phase there is a 
nonzero density of infinitely long defects (see figure 5b), 
which can be called a condensation of disclinations. 



VI. ORDER PARAMETERS 

An order parameter for a system with a local gauge 
symmetry must be gauge-invariant [29]. In our model the 
gauge group is Z2 and the local gauge transformation is 
given by eqn. (10). A product of spins and link variables 
is therefore gauge-invariant if the number of factors (Si 
or Uij) pertaining to each site i is even. Observables 
containing only gauge fields are therefore made of Wilson 
loops — products of the link variables around a closed 
curve 7 (typically a rectangle): 

W(j)= Y[ U u . (19) 
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(a) (b) 

FIG. 5. A caricature of the defect which is forced by bound- 
ary conditions to traverse the system in (a) the N or T phase 
and (b) the / phase. 

Purely spin-dependent quantities are gauge invariant if 
each spin enters an even number of times, such as the 
the defect counting term in eqn. (8), and the familiar 
traceless, symmetric, tensor order parameter 

Qa(3(i) = Si a Sip — S a p/3, (20) 

used in the Landau theory of the nematic-isotropic tran- 
sition. In the absence of an applied ordering field, long- 
range nematic order can be detected by calculating the 
correlation function 

N(i,j) = <Tr Q(i)Q(j)> = {(S; • Sj) 2 - 1/3). (21) 

N(i,j) is invariant under global spin rotation, and thus 
measures the tendency of the spins to align along a com- 
mon axis, regardless of its orientation. When Ql = 
lirri| r ._ r i^oq N(i, j) > 0, or equivalently (Q(i)) ^ 0, the 
system is nematically ordered. Thus, we have a local 
order parameter to distinguish the nematically ordered 
phase N from the isotropic phases / and T. 

Gauge invariant quantities can also be constructed by 
using both spin and link variables, as in the following 
path-dependent spin-spin correlation function: 

C(i, j; 7 ) = (s i -S J - JJ U kl ), (22) 

with the product of link variables taken along some path 
7 joining sites i and j. This is a generalization to arbi- 
trary gauge of the ordinary spin-spin correlation. In the 
infinite K limit, the choice of path is immaterial, since the 
product of links around any closed path is guaranteed to 
be one. Thus, in this limit, (C(i, j; 7)) is equal to (Sj-Sj) 
for the ordinary Heisenberg model at the same value of 
J. For any finite K, however, the path-dependent corre- 
lator (eqn. 22) always decays to zero exponentially with 
the separation between i and j, independent of the path 
7. This is not really surprising. For finite K, the links 
represent fluctuating degrees of freedom (the topological 



defects), and a single "weak link" will change the sign 
of the path-dependent product. (The same phenomenon 
occurs with the Wilson loop expectation value. Even at 
large K and in the presence of other ordering, it decays 
exponentially as a result of a non-zero density of very 
small defect loops.) 

An individual spin is not gauge invariant, so (S) 7^ 
is possible only if the gauge is fixed. Nevertheless, global 
rotation symmetry can be broken by a preferential align- 
ment of the spins along some axis. This is precisely what 
is measured by N(i,j) (eqn. 21) or (Q(i)) (eqn. 20). 
Ferromagnetic ordering implies nematic ordering but the 
converse is false. That the global SO (3) symmetry ac- 
tually is spontaneously broken at some finite, non-zero 
value of J for any K is clear from the results of the previ- 
ous section. We know that it is broken at finite J (points 
N and H on figure 3) in both of the limits K = 00 and 
K = 0. Reducing K introduces frustration and makes 
ordering more difficult. The only reasonable conclusion 
is that there is a line J C (K), monotonically decreasing 
with K, at which nematic ordering occurs. 

For the pure Ising gauge theory ( J = 0), the two phases 
on either side of the Ising transition discussed in section 
IV B above are distinguished by the asymptotic behavior 
of large Wilson loops (W(L,T)), for which the path 7 is 
a closed rectangular loop of sides L and T. In the small 
K phase, the Wilson loop follows an "area law" for suffi- 
ciently large L and T: \n(W(L, T)) is proportional to the 
area LT of the closed path. This corresponds to the "con- 
fining" phase of the gauge theory. In the large- K phase, 
on the other hand, sufficiently large Wilson loops obey 
a "perimeter law," so that \n(W(L,T)) is proportional 
to the perimeter L + T of the closed path. This corre- 
sponds to the "free-charge" phase of the gauge theory. 
These results emerge from expansion methods described 
further in section VII A. As also shown there, the "defect 
line tension" ST/L is related by duality to the spin-spin 
correlation function of the Ising spin model. 

The distinction between the fully disordered and topo- 
logical^ ordered isotropic phases is more subtle when J 
is nonzero. There is no local order parameter which dis- 
tinguishes between the I and T phases. The line tension 
is therefore quite important because it does perform this 
function. Coupling of spin and gauge degrees of freedom 
results in a decay of Wilson loops which is asymptoti- 
cally a perimeter law for any nonzero J, with a crossover 
scale of L x = 41og(J/3)/log(tanh K). The defect line 
tension described in section V, however, remains a valid 
diagnostic of topological order. 

The results discussed so far show quite clearly that 
there is a line connecting the first-order transition at N 
to the second-order one at H as depicted in figure 3, at 
which global spin rotational symmetry is broken, there 
is also a critical point at I. We have until now somewhat 
implicitly assumed that this point actually separates two 
completely distinct phases / and T, though the line ten- 
sion appears to be a good diagnostic of this distinction, 
and offers strong support for it. Accepting this point 
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(which is established more carefully in later sections), the 
simplest possible phase diagram topology, i.e. that with 
the fewest number of phases, is shown in figure 3, which 
in fact is the outcome of our Monte Carlo simulations. 

In section VII, we demonstrate analytically that the 
continuous Ising transition near I persists for nonzero J, 
and the Heisenberg transition near H persists for finite 
defect activity e~ 2K . In section IX we fill in the ana- 
lytically intractable interior of the phase diagram using 
Monte Carlo simulation. This work shows that the three 
transition lines emanating from the border of the phase 
diagram meet at a multicritical point M. The jump of 
the order parameter at the symmetry-breaking transition 
goes to zero at M, so that the transition is first order to 
the left (smaller K) of M and continuous to the right. 
Finite size scaling has been employed to verify that the 
entire continuous transition line is in the Heisenberg uni- 
versality class, not just the point at infinite K. Calcula- 
tions of the specific heat strongly suggest the existence of 
the Ising transition line originating from the pure gauge 
theory transition. This is verified by calculations of the 
defect line tension. 



VII. ROBUSTNESS OF CONTINUOUS 
TRANSITIONS 

In this section, we demonstrate that the continuous 
Heisenberg N/ T and Ising T/I transitions found in sec- 
tion IV at K = oo and J = respectively, persist for 
K < co and J > 0. These transitions are therefore 
generic and occur for a finite range of material parame- 
ters, and not just at isolated points. Our tool is pertur- 
bation theory in e~ K and J, which can be developed in 
terms of polymer expansions. The general formalism of 
such expansions is developed in subsection VII A. Sub- 
sections VII B and VII C treat the limit K — ► co, while 
subsections VII D and VII E deal with the limit J — > 0. 
Finally, we treat the limit K — ► in subsection VII F. 



factor arising from the spin term of our Hamiltonian, eqn. 
(11). Expansion of the ensuing product produces a sum 
of terms, each of which contains factors Uij tanh(JS;-Sj) 
for the links (ij) in a distinct set. We decompose each 
such set into constituent "polymers" connected at the 
lattice sites, and define the activity of a polymer co as 



p{uj) 



(24) 



where dSi denotes the usual integration measure over the 
directions of S; (normalized: J dS = 1), and 



p(u,S) = Yl tanh(JSjS fc ) 



(25) 



(jk)eu 



is a quantity we will need later when discussing small K . 
A polymer clearly has zero activity unless it is closed, 
i.e. an even number of constituent links impinge on each 
lattice site. (When calculating a correlation function in- 
stead of the partition function, we wil alter this definition 
slightly, so that some open polymers may have nonzero 
activity.) In this representation, the complete partition 
function is written 



Z(J,K) = Y J Z{K)(W{C)) G J] pM, 



(26) 



u>EC 



where the sum runs over collections C of non-intersecting 
spin polymers, and W(C) is the associated generalized 
Wilson loop, a product of factors Uij for each link in C. 
Its expectation, evaluated here in the pure gauge model 
at coupling K (as indicated by ()g), can also be written 
as 



(W(C)) C 



Z C {K) 
Z(K) ■ 



(27) 



A. Polymer Expansions 



where Z(K) is the gauge theory partition function at 
coupling K, and Zc(K) the partition function in the 
presence of "source" C. 

Spin polymers interact via both the gauge field and a 
hard core repulsion. 



Polymer expansions for lattice models express the par- 
tition function and correlation functions in terms of an 
interacting system of self-avoiding chains of links on 
the lattice. The simplest example is the familiar high- 
temperature expansion of the Ising model; Mayer expan- 
sions may also be viewed as polymer expansions. For 
more details on this subject, the derivation of eqn. (30), 
and careful discussions of the convergence question, the 
reader is directed to references [30-32] . 

With Ising variables (a = ±1), the simple identity 

e xa = cosh(x)(l + a tanh(x)), (23) 

is remarkably useful. By taking a = U t j and x = JS t ■ Sj, 
this identity can be applied to the part of the Boltzmann 



B. Large K 

Now we evaluate (W(C))g for some fixed C, by intro- 
ducing a second type of polymer. At large K, the gauge 
field configurations are most conveniently represented in 
terms of defect loops on the dual lattice. A link on the 
dual lattice pierces a unique plaquette on the original lat- 
tice. For a configuration of {Uij}, that (dual) link is part 
of the defect network if the product of U^s around its 
associated plaquette is —1. As with the spin polymers, 
we decompose the defect network into pieces which are 
connected at the dual-lattice sites. We call the pieces 
"defect loops," although the nomenclature is not ideal 
since such a loop may cross itself many times. 
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A defect loop 7 of total length |7| has a C-dependent 
activity given by 

z( 7) C) = (_ 1 )i(7,c) e -2K| 7 | ! (2g) 

where 4(7, C) = ±1 is the parity of the linking of 7 with 
C, which is equal to the product of the linking parities of 
7 with the separate constituents of C. (For two loops 7 
and u>, i(j,uj) is —1 if 7 wraps around u> an odd number 
of times, and +1 if an even number). 

By construction, the defect loops do not overlap. We 
may also think of this as arising from a hard-core repul- 
sion. The partition function in the presence of source C 
will then be written as 

w=En^ c ) n a-'fr.v)), (29) 

D 7G-D 7,7'GD 

where £(7,7') is +1 if the defects 7 and 7' overlap, and 
is if they don't. The sum over collections D of defects 
does not then need to be restricted to non-overlapping 
sets. Each term in the expansion of the product of fac- 
tors (1 — Z(7, 7')) can be associated to a graph whose 
nodes represent the polymers 7, and in which two nodes 
corresponding to 7 and 7' are joined by a line if they 
overlap, i.e. if Z( 7 , 7 ') = 1. 

This formulation is advantageous when we pass to the 
logarithm, i.e. the free energy. Explicitly, 

F C (K) = -\nZ c (K) = J^J R zfr), (30) 

D ' ' ' -yED 

where the sum is over connected graphs D containing \D\ 
nodes. The "index" of D is defined as 

n(D) = ^'(-l)'(G), (31) 

GcD 

where the sum is over connected subgraphs G C D which 
contain all the nodes of D (thus \G\ = \D\), and 1(G) is 
the number of lines in G. The formula (30) for the free 
energy is an expansion (cf. eqn. 28) in powers of 

a ee e- 2K . (32) 

For a sufficiently small (K sufficiently large), the com- 
bined influence of factors of a and the requirement of 
connectedness suppresses the higher-order contributions 
enough that we can prove convergence of the expansion. 

We can now see that \n(W(C)) G = In Z c (K)-\n Z(K) 
is a difference of two expansions like eqn. (30), one with 
activities z(-y,C) and the other with z(j). Only terms 
containing polymers with 21(7, C) ^ 21(7) need to be cal- 
culated; all others are killed by the subtraction. To see 
this in action, we calculate the contribution of the small- 
est defect loops to (W(C))a- In that case, the D's oc- 
curring in the sum in eqn. (30) consist of single defect 
loops of four links going around the perimeters of elemen- 
tary plaquettes on the dual lattice. Only if 7 encircles C 



are z(-y,C) = —2(7) 2(7) = a 4 different. The result is 
therefore 

\n(W(C)) G = -2e- 8K ■ \C\ + 0( e - 12K ). (33) 

This result can now be used to derive an effective 
Hamiltonian for the spins alone. We sum over the de- 
fect configurations to express Z(J,K) purely in terms 
of spin-polymers. In the expansion of eqn. (26), each 
factor tanh(JS; • Sj) in a spin-polymer activity p (eqn. 
24) is accompanied by a factor of Uij, which goes into 
W(C). From eqn. (33), we can think of each link in C as 
bringing a factor exp(— 2a 4 ) to the expectation (W(C)). 
Alternatively, we can put this factor with the tanh to 
write, correct to order a 4 , 

Z{J,K)= /"jJdSi JJcosMJSj-S*) 

i (jk) 

x (l+ er 2 " 4 tanh( JSj ■ S*)) • (34) 

As remarked before eqn. (24), the constraint to closed 
loops is still in force by virtue of the integration over the 
spins. 

Taking a logarithm gives an effective Hamiltonian 
H * s = J]{(l-2a 4 )JS i -S i 

(ij) 1 

+ a 4 f] ( " 2J8 ;,' Sj ' )n } + 0(a«). (35) 

n=2 ' ' 

To this order, the effective Heisenberg coupling is reduced 
to J e ff = (1 — 2a 4 ) J. Since eqn. (35) contains additional 
terms higher order in the spins but the same order in 
a, we can say J c (a) — J c (a = 0) = (D(a 4 ), but cannot 
predict the numerical prefactor. 

The effective Hamiltonian in eqn. (35) gives the same 
free energy as the original model. If one wants to com- 
pute expectations, it is helpful first to restore gauge in- 
variance by multiplying S t ■ Sj by Uij, where these link 
variables are subject to an infinite effective K. Going 
through the same arguments, one finds that expectations 
of quantities involving link variables Uij (the path depen- 
dent scalar product C(i, j; 7 ) of eqn. 22, for instance ) 
can be computed by replacing Ui 3 inside the expectation 
with e~ 2a U^ and then using H e g . Thus, exponential de- 
cay of C(i,j; 7) is maintained even in the ordered phase. 

The program of eliminating the link variables U ap- 
pears to be going well, and the result (eqn. 35) argues 
for suggestive of the irrelevance of small defect activity. 
However, this is the lowest order result; only the smallest 
defects have been eliminated. Higher order (in a) terms, 
while individually small and formally irrelevant, prolifer- 
ate alarmingly at higher order. It is not entirely obvious 
that they can be neglected, or that the "sum" even ex- 
ists. In section VII C below, the issue of defect activity 
irrelevance is taken up again, within a renormalization 
group approach. 



11 



r iW 



•— e- 



1 



• • J« — • — '<> 



• • ii • it 



• • ] 



(a) (b) 
FIG. 6. The lattice paths which occur in the Real-Space 
blocking scheme for (a) the spins and (b) gauge (link) fields, 
in the case of two dimensions. 

C. Irrelevance of Defect Activity 

The first task in developing a real-space renormaliza- 
tion group approach is to define block variables. We need 
to block both spins and link variables. The spins can be 
handled in a usual way, with a bit of care to insure gauge 
invariance. The link variables are not so familiar in this 
context and require some thought. We will present one 
particular blocking scheme, and then comment on the 
motivation. 

Figure 6 illustrates the blocking procedure. Sites in 
the blocked lattice are denoted with a prime superscript, 
thus, site i' is associated with the cubical block B(i') of 
sites centered on i' with side length 2L + 1: 

fl(i') = {j : \x a (j) ~ x a (i')\ < L, a = 1, 2, 3}. (36) 

The contribution of the spin at j to the block spin S[, 
is found by parallel transporting it to i', using some set 
of standard paths. A convenient choice, illustrated in 
fig. 6 for the case of two dimensions, is to move along 
coordinate directions in reverse lexicographic order, i.e. 
move along the z-direction until the z-coordinate equals 
that of the destination, then along the y-direction, then 
along the x-direction. Denote this path from j to i' by 



The blocked spin is now defined as 



S' 



where the parallel transport factors are 



t/(r^)EE ][ u kl . 



(37) 



(38) 



The rescaling factor (s will be adjusted to keep the dis- 
tribution centered around vectors of unit length. 



For two neighboring sites i' and j' in the blocked lat- 
tice, we might try to make a blocked link variable by sim- 
ply taking the product of link variables along the straight 
path from i' to j' denoted by fV j< . This is a really a dec- 
imation scheme. It turns out to be much better to take 
products of link variables along paths which are defor- 
mations of this straight path on a scale L and average 
them. Defining closed loops 

^i',k,i,f = Fi',k U Vk,i U Vj^i U r^j/, (39) 

our blocked link variable is (see figure 6) 

ul, j ,=CuU(t i .j.)jjY l 'u(r i ',k,ij)- (4°) 



The prime on the sum denotes restriction to k G B(i'), 
I € B(j'). Again, there is a rescaling factor to keep the 
link-magnitude distribution centered at one. 

The preservation of gauge invariance by our blocking 
scheme is easily checked. Performing the gauge transfor- 
mation Si — ► 4>iSi, Uij — ► 4>i(j)jUij on the original lattice 
first (<f>i = ±1 chosen independently for each i), and then 
blocking results in 



s 4>i- s 



u', 



(41) 



This is a gauge transformation on the block lattice, as 
required. 

After m renormalization group steps, the running 
Hamiltonian looks like 



ff m = £* m (S?) + £*£(tf?.) 

+ Jm E UijSi ■ Sj + K m JJ UijUjkUkiUu + 



(42) 



and we need to determine the flow of the various cou- 
plings. The first two terms provide the weighting for 
the magnitudes of the spins and link variables and their 
indicated functional dependence is dictated by gauge in- 
variance and rotational symmetry. On the original lattice 
(to = 0), they were delta functions at one. The ellipsis 
indicates all other interactions generated by the renor- 
malization group transformation. 

Our blocking scheme treats large and small defect loops 
differently, and this is reflected in the structure of H m . 
Defects smaller than the current lattice spacing L m are 
effectively integrated out in the sense that they no longer 
appear as extended objects, but they suppress the mag- 
nitudes of the blocked links. Their effects are reflected 
in After blocking, but before rescaling the link 



variables, Uf,-, 
1 j 



typically has a value ~ e 8£ ", from the 



small defects, (a = e~ 2K ) This has an influence of the 
magnitude of the product Yl^' j'£> 'j'k'Uk'i'Uini around a 
blocked plaquette. That product, like the blocked link 
variables Uv making it up, is an average over many 
paths. It will not be negative unless most of them are 
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encircled by a defect. This can be accomplished by one 
large defect of length ~ L, or an even greater length of 
small ones. 

The value of K m is therefore determined by the proba- 
bility to have such a large length of defect (~ L m in orig- 
inal lattice units) through the blocked plaquette; every- 
thing else is absorbed into As in the calculation of 
A'eff in section VII D, the spins also contribute to the 
renormalization of K, and tend to suppress it. Let us 
neglect this, though, and make the conservative estimate 

K m+1 ~ L(K m -c)e- 16La ™, (43) 

where the c comes from an entropy factor (the number 
of distinct defects of total length L is ~ expcL). The 
rescaling of the link variables has also been taken into 
account in making this estimate and accounts for the 
final exponential factor. 

As long as the initial value of K is large enough, and 
we don't try to take too large blocking steps, the defect 
fugacity a is being rapidly driven to zero, exponentially 
fast, in fact. If K starts out too small, the entropy dom- 
inates, and this scenario no longer holds [33]. 

Thus, in order to conclude that a small defect activity 
does not alter the Heisenberg universality of the transi- 
tion, it is only necessary to verify that in the process of 
vanishing it does not induce any extra interactions which 
are themselves relevant. The fundamental constraints 
on the form of these are gauge and rotational invari- 
ance. Among terms involving only the gauge variables, 
there are even more irrelevant multi-plaquette terms, 
and also the bond weight § u , which ultimately amounts 
to annealed randomness in the Heisenberg model bond 
strengths. Two of the least irrelevant pure-spin terms 
which can occur are (Si • Sj) 2 and the spins-around-a- 
plaquette term Yl(St ■ Sj), which will be seen again in 
section VII F. None of these perturbations is relevant. 
Thus, the presence of weak defect activity does not alter 
the (Heisenberg) universality class of the transition. 

D. Small J 

At small J, the spins are a small perturbation to the 
pure gauge theory. We can try to integrate them out, 
as we did the gauge field excitations in section VII B, to 
obtain an effective Hamiltonian for the gauge field, with 
corrections ordered as a power series in J. The polymer 
expansion is (see 24,26) 

z(j,k)= ^2 ^ flG En (pm n u *j) • ( 44 ) 

{Uij} C ujGC \ (ij)doj J 

At leading order, the only spin polymers occurring are 
those which go around an elementary plaquette on the 
lattice, just as with the defect loops on the dual lattice 
in the K — ► oo limit. The required activity is 



p= [dS] J] tanh(JSi • Sj) 

= 3(J/3) 4 + 0(J 6 ), (45) 

where the product is over the links around the plaquette 
P. Exponentiating this result gives an effective action 
for the gauge field [34] with A' eff = K + J 4 /27 + £>( J 6 ). 

Terms of higher order in J also introduces gauge cou- 
plings on larger closed loops of links which decrease ex- 
ponentially with the length of loop. 

The polymer expansion for the spin problem is known 
to be convergent for small enough J, so the only question 
is whether it can legitimately be rewritten as an effective 
action for the gauge field. The entire calculation is very 
similar to the one performed for the large K regime in 
section VII B, with the defect loops replaced by spin- 
polymers. In fact, if our spins were Ising spins, there 
would be a perfect duality. We believe that the near- 
duality makes this effective action calculation as solid as 
the previous one, and establish the robustness of the Ising 
transition for J > 0. 



E. Line Tension Again 

The defect-line tension was discussed heuristically in 
section V. Additional insight can be gained by using 
the duality between the three-dimensional Ising spin and 
gauge models. This duality amounts to no more than the 
banal observation that the polymer representations of the 
ordinary nearest-neighbor Ising model H = —J^2<Jia 3 
and the gauge model in eqn. 11 are identical (up to an in- 
nocuous overall factor of (cosh J) 3N ) when the couplings 
are related by 

e - 2K = tanh J. (46) 

The self-duality of the Z 2 gauge- Higgs model (i.e. 
spins and links are Ising variables) is just as clear. Re- 
ferring to section VII B, we see that the spin polymers 
and defect loops have identical activities apart from a 
linking-number related interaction, which is completely 
symmetric between the two types. 

Dual representations of correlation functions are also 
easily written down by going through the polymer repre- 
sentation. For instance, (<Ji<Jj) -Z corresponds to polymer 
configurations having one polymer with an odd number 
of links connected to sites i and j. One way to generate 
this set is to take one fixed configuration S which satisfies 
the condition, a single path of links between the two sites 
for instance, and then forming the symmetric difference 
with each configuration for Z, which are made of closed 
loops. (The symmetric difference of two sets A and B is 
AAB = (AuB) \ (A n B).) 

In the gauge model, the set of links dual to plaque- 
ttes on which the field strength is —1 is required to be 
closed by the nature of the underlying fields. If the cou- 
plings are not all the same, the set of plaquettes on which 
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KpUp, i.e. on which the configuration is an excitation 
over the local ground state is not required to be closed. 
The defect loops must really be identified with the ex- 
cited plaquettes, so if we reverse the sign of the coupling 
on plaquettes in £*, with £ as above, the partition func- 
tion of the resulting model is the dual representation of 
Z {<Ji<Jj). We get exactly the same thing by evaluating 
the expectation value of the operator 

£>(£*) = Y[ exp(-2K Up) (47) 

Pes* 

with the original Hamiltonian. This disorder operator, 
(a 't Hooft operator [35] for our particular gauge group) 
is therefore the dual representation of the operator <Ji<jj. 
The number of defect lines entering or leaving any ele- 
mentary cube of the lattice has the same parity as the 
number of plaquettes on that cube for which the coupling 
has been reversed in sign. 

The 't Hooft disorder operator can detect the tran- 
sition in the pure gauge theory (J = 0). As i and j 
become infinitely separated, (£)(£*)) tends to a non-zero 
constant at small K and decays exponentially in \i — j\ 
for large K. At J = 0, the only relevant aspect of the 
path £ is its endpoints. For J > 0, this is no longer 
true, since the spins are sensitive to the signs of Up and 
not just the products KpUp. As a result, for J > 0, 
(£>(£*)) decays exponentially in the length of the path 
|£| for any value of K. This behavior is reminiscent of 
the path-dependent spin-spin correlation function (eqn. 
22). This is not an accident; in the case of 7Li spins, the 
two are dual to each other. 

The defect line tension is similar to the 't Hooft op- 
erator, and equivalent at J = 0, where the the tension 
is therefore positive for K > K c and zero for K < K c . 
Unlike the Wilson loop or 't Hooft operator, however, the 
line tension continues to be a good diagnostic for nonzero 
J. Rigorous results [36] show that it vanishes in some re- 
gion around J = K = and is strictly positive in some 
region around K = oo, J = in our phase diagram. 
Whether these behaviors hold throughout the entire I 
and T phases however, is not rigorously established. 

To go further, consider the truncated energy-energy 
correlation 

(e(ij); e(kl)) = (ai(Tj(Ti(T m ) ~ (<Ti<Tj)(<Ti<r m ), (48) 

for nearest neighbor pairs (i, j) and (l,m), which yields 
the specific heat when summed over links (Im). Choose 
£ (as in the spin-spin correlation construction above) to 
be the two plaquettes dual to the links (ij) and (Im), so 
that the dual representation of (e(ij) e(kl)) is obtained 
by weighting each configuration of defects contributing 
to the gauge theory partition function by an extra factor 
of e~ 2K for each of those links which does not occur in the 
defect set and e +2K for each which does, and similarly 
for (e(ij))(e(kl)). 

We use the notation x(u) = U^y, f° r the product of 
link variables around the plaquette (ij)* which is dual to 



the link (ij). Thus, x(ij) = 1 if dual link (ij) is in the 
defect set and x(*j) = otherwise. Then we have the 
correspondence 

(e(ij); e(kl)) 
= 4sinh 2 2K( X (ij);x(kl)) 
= smh 2 2K(U itj) ,;U {kl) «} (49) 

The divergence of the specific heat in the spin model 
shows that this function becomes long-ranged at the crit- 
ical point. In the gauge theory language, what eqn. (49) 
measures is the degree to which frustrations of widely 
separated plaquettes are not independent. We argue that 
its asymptotic behavior reflects that of the probability 
for two widely separated points to be in the same defect 
cluster. Multiplying eqn. (49) through by a factor of the 
partition function Z, we rewrite the result as a sum over 
the defect cluster [37], 7, which contains dual link (ij), 
since only such configurations contribute. We will also 
extract an explicit factor 21(7) = ck'" 4 ' of the weight of 
cluster 7 (eqn. 26 with C = 0). Thus we write 

£'z( 7 ) (1 - (x(kl)) A ) + «x(*0>AVr " <x(*0>a) , 

7 7 

(50) 

where the first sum is over defect clusters 7 which contain 
both (ij) and (kl), and the second over those containing 
(ij) but not (kl). The subscripts on the expectations 
indicate the lattice we are calculating for - the full lattice 
A, or with the cluster 7 removed, A \ 7. 

If K is small enough, the probabilities of (kl) G 7 and 
of |7| > M asymptotically decay exponentially with the 
distance between (kl) and (ij) or M , respectively. Sim- 
ilarly, the difference (x(kl)) — (x(^O)a decays exponen- 
tially with the distance of (kl) from the boundary of A 
(the former expectation is for an infinite lattice with no 
holes). If the correlation eqn. (50) is to exhibit power law 
decay at K = K c , at least one of these three quantities 
must also show such a change in asymptotic behavior. 
It seems extremely unlikely that the effect of cutting out 
small pieces of the lattice has a qualitatively slower falloff 
than the cluster size. Accepting that, the transition at 
K = K c , J = is accompanied by a percolation of defect 
lines. The difference from ordinary bond percolation is 
that the frustration network is made of closed loops and 
thus is not allowed to have free ends. 

This picture is easily related to the line tension. The 
correlation function of eqn. (49) and our manipulation of 
boundary conditions in section V both serve to to mea- 
sure the ease with which a defect line can join two distant 
points. The line tension is the inverse of the correlation 
length for the frustration-percolation. 

This percolation criterion must hold for J > to be 
a useful description of the // T transition. There is little 
doubt that it is more robust against the perturbation 
of positive J than either the Wilson loop or 't Hooft 
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disorder operator. A crude estimate says that for e~ 2K < 
J, the leading behavior of the large- A defect expansion 
is the same as at J = 0. For the small- A expansion 
in terms of plaquette-surfaces, the leading terms involve 
a tube running between the two plaquettes in question 
and the J perturbation is negligible for tanh A > J 2 . 
As stated in section V, the line tension cannot vanish 
in phase N, since that would certainly destroy the long- 
range ordering. The case of the Z2 Higgs-gauge model 
[38], however, argues for some caution. It has only two 
phases - a free-charge phase at large A and small J, and 
a confinement-screening phase everywhere else. The line 
tension, which vanishes in some region near K = J = 0, 
becomes positive at large J without passage through a 
bulk phase transition. 



F. Small K 

In section IV, we showed that for K = our gauge 
theory becomes a spin-only model akin to the Lebwohl- 
Lasher model. Now we will use the methods introduced 
in the previous section to extend the elimination of the 
gauge variables to small non-zero K . 

We expand the spin part of the Boltzmann factor as 
before, but without integrating over the spins, and use 
the pure gauge Hamiltonian to evaluate the expanded 
form term- by-term. Explicitly, 

Z = j[dS}Z[S], (51) 

with 

Z[S] = JJcosh(JSi • Sj) 

(ij) 

x / JJ(1 + Uij tanh( JS, • Sj))\ . (52) 

Equivalently, this can be written as 

Z[S] = Z[S;K = 0]J2Zc(K)l[p(u J ,S), (53) 

C ujGC 

where the sum is over collections C of closed graphs on 
the lattice, as in equation 26 for the large A case, and 
p(u>,S) is defined in eqn. (25). The factor Zc{K) is 
the pure gauge theory partition function in the presence 
of source C (eqn. 27). The gauge field induces a weak 
non-contact interaction between the spin polymers. 

By use of eqn. (23), the Boltzmann factors associated 
with the pure gauge theory can also be rewritten as 

e KUp = cosh K(l + Up tanh A"), (54) 

where we use the shorthand Up = Yi(ij) Uij for the prod- 
uct of link variables around a plaquette P. Expanding 
the product gives us collections of plaquettes with fac- 
tors of tanh A'. The expansion of the spin part is as 



before. Upon summing over ±1 for each link variable 
Uij, any surviving term of an expansion must have an 
even number of occurrences of each link variable. Thus, 
our polymers will consist of edge-connected collections of 
plaquettes each carrying a factor tanh A', and a factor 
tanh(JS;Sj) for each boundary link (ij), if any. In the 
presence of a Wilson loop, the loop and all the plaquettes 
from expansion of the gauge Hamiltonian or links from 
the spin Hamiltonian which overlap it on a link are to 
be considered as a single polymer, whose activity is zero 
unless all the factors of Uij are cancelled in pairs, and is 
otherwise evaluated as for the others. 

The leading terms contributing to a Wilson loop in 
this case are (i) the one containing all the plaquettes on 
a surface bounded by the loop, which gives an area law 
decay in the case of J = 0, or (ii) the spin polymer which 
tracks along the Wilson loop, which gives a perimeter 
law and is dominant for large enough loops at J > 0. In 
either case, what determines the minimal polymer is the 
need to cancel the factors of U along the loop. In the 
expansion of eqn. (53), the expectations are in the pure 
gauge Hamiltonian, so there is area law decay and the 
expansion is well under control. As a first approximation, 
we keep only single plaquette graphs. Exponentiating, 
the result is 

^iff =E 1 °g[ C0Sh ( JS ^ S J )] 

+ J'J2 II tanh ( JS fc -Si). (55) 
□ fczen 

Here, J' = J 4 tanh K. Notice that the correction term 
in equation (55) is almost the same as the disclination 
counter equation (8), just without the sign function. 

VIII. OTHER PERTURBATIONS 

The model we have been considering contains some ex- 
act symmetries which are only approximate for real ne- 
matics. These are (1) local head-to-tail symmetry, i.e. 
S — ► — S, and (2) global space-spin rotation invariance, 
under which all the spins are subjected to the same ro- 
tation ( S(r) — y RS(r)). Local head-to-tail symmetry is 
lacking despite the fact that global head-to-tail symmetry 
is unbroken, i.e. that the heads and tails do not order. 
The absence of a globally broken symmetry does not im- 
ply that the Hamiltonian lacks terms breaking it locally 
- only that they are too weak to induce long range order 
(vector ordering in this case). 

The absence of global spin-rotation invariance is ex- 
hibited explicitly in the Frank free energy (eqn. 1) for 
generic A; (it is invariant at the special point Ai = A 2 = 
A 3 ). Since the spins of our model are related to the phys- 
ical orientation of extended bodies, it is not surprising 
that the only exact symmetry is that under simultane- 
ous and identical rotations in spin- and real-space. 
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We show here that small perturbations of the Hamil- 
tonian which break these symmetries do not change our 
results, so that we can feel safe applying them to real- 
world nematics. 



A. Lack of Local Inversion Symmetry 

First, we consider the local head-to-tail symmetry. We 
need an extra variable on = ±1 for each site which keeps 
track of whether S t is oriented with the local vector. It 
has identical gauge transformation properties to that of 
Si itself, so that the product ot{Si is the "real" local vec- 
tor (and gauge invariant as it must be). A direct Heisen- 
berg interaction is added to the Hamiltonian of eqn. (11) 
by the term g Yl(ij) onajSi-Sj, for some small g (compare 
the "Mattis spin-glass"), where the a's are to be summed 
over in the partition function. As far as the spins S are 
concerned, this is equivalent to adding a second indepen- 
dent gauge field with K = oo, thus the Heisenberg char- 
acter of the transition is clearly preserved. It is also easy 
to see that the interaction does not produce a true vector 
ordering once the nematic phase is entered. The tendency 
toward complete vector ordering is only strengthened if 
all the spins are forced to orient along a common axis, 
yielding an ordinary Ising model in the variables c^S^. 
But the coupling g is weak, so there will be no ordering. 

B. Space-spin Coupling 

Near the N/ Ttransition we have seen that our model 
can be mapped onto an effective Heisenberg model. In 
the critical regime, a familiar transformation [1] converts 
this fixed spin model into an n = 3 "soft-spin" |S| 4 Lan- 
dau theory with Hamiltonian density 

^ = ||VS| 2 + | r | S | 2 +^|S| 4 , (56) 

where S is now a three- vector with unconstrained length. 
To incorporate space-spin couplings, we introduce terms 
which mix up the spatial component indices and the spin 
components. The only such terms with any chance of 
altering the critical behavior have as few gradients and 
powers of S as possible. Indeed any term involving more 
than two powers of each will be strongly irrelevant. The 
only potentially relevant term is therefore |V • S| 2 . This 
perturbation was analyzed [39] in the heyday of RG in 
an e = 4 — d expansion, and found to be irrelevant, with 
renormalization group eigenvalue A = — e 2 /108 + 0(e 3 ). 
Thus, there is nothing to fear from space-spin couplings, 
either: the Heisenberg universality class of the transition 
survives. (This issue is discussed further in Paper II.) 



IX. MONTE CARLO RESULTS 
A. General Methods 

We have employed Monte Carlo simulation to investi- 
gate the lattice gauge model of eqn. (11). The vast ma- 
jority of runs were for three-component spins, though we 
have also investigated four-component spins. The sim- 
ulations were implemented on a Sun SparcStation with 
a standard Metropolis algorithm on lattices of size up 
to 16x16x16 sites. Most of the runs employed periodic 
boundary conditions for all variables (spin and gauge) 
to eliminate boundary effects. A combination of free 
and fixed boundary conditions, however, was necessary 
to measure the line tension (see section IXC and V). 

Instead of allowing the spin variables to sample the 
entire unit sphere, we used a discrete set of allowed val- 
ues. This simplifies and speeds the simulations, since 
Boltzmann factors can be stored in a lookup table, and 
it is also somewhat easier to select candidate Monte 
Carlo moves. For the three-dimensional order parameter 
space, we used the most symmetric set of allowed vec- 
tors, namely the 30 vectors pointing to the centers of the 
edges of an icosahedron. Such a discretization represents 
an anisotropy of the single-spin weight which breaks the 
original 0(3) rotation invariance to its largest discrete 
subgroup, the icosahedral group Y. 

In order for this anisotropy to be irrelevant at the 
Heisenberg critical point (in the renormalization group 
sense), the allowed spin orientations must to cover the 
sphere sufficiently uniformly. The icosahedral edge vec- 
tors easily pass [40] this test. Use of a coarser discretiza- 
tion which is still irrelevant presents two problems. First, 
it will be necessary to get closer to the transition be- 
fore crossing over to the fully-symmetric critical behav- 
ior. Secondly, the discretization introduces a spurious 
freezing transition, due to the presence of a spin wave 
gap. It is desirable to push this artifact deep into the ne- 
matically ordered phase. The same effect also results in a 
small shift of the numerical value of the critical coupling. 

Among the measurables which were extracted from the 
simulations were 

1. the average plaquette value P = (UijUjkUkiUu), 

2. the gauge-invariant nearest-neighbor correlation 
C = (UijSi ■ Sj), 

3. average Wilson loop of several sizes, and 

4. a scalar measure of the strength of nematic ordering 
given by 

2 N /3 m 1 \ 

where Q a p = (l/N)^2 i Q(i) a p is the nematic order pa- 
rameter, and N is the number of sites in the lattice. 
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Tr Q 2 is invariant under rotation and measures the de- 
gree of alignment of the individual spin axes for a uni- 
axially ordered phase. The normalization of eq. (57) is 
chosen so that q 2 vanishes in the disordered phase and 
is equal to unity in the completely ordered state at zero 
temperature. Terms with i = j are excluded from the de- 
finition of q 2 ; they would make a constant contribution 
subleading in 1/N. Leaving them out ensures q 2 = in a 
fully disordered state. Normalization is such that q 2 = 1 
in a completely ordered state. 

We studied the equilibrium values of these quantities 
by stepping along lines in the J — K phase diagram at 
a variety of orientations, allowing the system to reach 
equilibrium after each incremental change, and then mak- 
ing measurements. The "temporal" development of these 
quantities for values of the couplings close to the transi- 
tions were examined by eye to determine the equilibra- 
tion time. Thermalization typically required 700-2,000 
update attempts per degree of freedom (link or site) with 
each coupling step. Typically 6,000 to 15,000 indepen- 
dent measurements were made at each coupling value, 
spaced by two or three update sweeps through the lat- 
tice to give reasonable statistical independence. This was 
sufficient to extract equilibrium averages of the quantities 
(1-4). To obtain accurate results on fluctuations about 
equilibrium values (e.g. the specific heat), as well as to 
examine finite-size scaling, 30,000-60,000 measurements 
were required at each coupling step. 

Our data confirm the existence of the three phases 
which were predicted analytically in previous sections. 
They also show that the transition between phases A^and 
/is first-order, and the transitions between A^and Tand 
between / and Tare continuous. As is well known, how- 
ever, it can be difficult to definitively establish the order 
of a transition via Monte Carlo. We have used finite-size 
scaling for the continuous transitions and shown phase 
coexistence at the first order N/I boundary. 

The distribution of Tr Q 2 for values of J near the 
N/I transition shows it to be first-order at K = 0. 
The double-peaked distribution shown in figure 7 demon- 
strates the coexistence of ordered and disordered phases. 
Clear evidence of a discontinuity in the order parameter 
expectation value q 2 is also seen, a feature which becomes 
sharper with increasing lattice size. 

Figure 8a displays the development of the order pa- 
rameter q 2 near the ordering transition close to the 
point where it meets the Ising transition line. The data 
strongly suggests that the nematic ordering transition is 
first order to the left of that point and continuous to its 
right. The mean plaquette value and specific heat also 
show qualitatively different behavior from one side of that 
point to the other. The difference is discernible even be- 
tween K = 0.73 (on the first order side) and K = 0.77 
(continuous). On the first-order side, the specific heat 
peak sharpens rapidly as the lattice size is increased, 
compared to a much more gradual sharpening for larger 
K at which the transition is continuous. This behavior of 
the specific heat is shown in figure 8b (see also the next 
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FIG. 7. Probability distributions for order parameter at 
K = in a 12 3 lattice for several values of J near the 
first-order transition. As J increases, the development of a 
bump corresponding to the ordered state and consequent di- 
minishing of the disordered-state bump at zero order para- 
meter is evident. In the thermodynamic limit, at any given 
J, only one of the bumps is expected to survive. The distri- 
butions were constructed by making histograms of a Monte 
Carlo "time" series. 
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FIG. 8. Development of the nematic order parameter 
Tr Q 2 at the transition, and near the crossover point from 
first-order to continuous. All data are for a 12 3 lattice. The 
apparent sharpness of the jump for K = .73 is essentially ac- 
cidental; neighboring data points straddle the region where 
Tr Q 2 begins to rise from zero. 

section and figure 10). 

B. Finite-Size Scaling 

Finite-size scaling analysis is a well-developed tool for 
the determination of critical exponents and orders of 
transitions (for reviews, see references [41] and [42]). The 
finite-size scaling ansatz for critical points asserts that in 
finite geometries characterized by a linear dimension L, 
scaling forms depend on a scaling variable L/£, in addi- 
tion to those already present at infinite system size. For a 
magnetic system, the critical coupling can be determined 
from the ratio 

ML,t)= { j^*f n (Ln (58) 

Here, m is the spontaneous magnetization, and n is a 
small integer. The indicated functional dependence fol- 
lows since both numerator and denominator have the 
same scaling dimension. Thus, at the critical temper- 
ature, this quantity is independent of L. Furthermore, 

lhm|/(LOocL^, (59) 

as can be seen by noting that this derivative must be fi- 
nite and non-zero as t — ► at fixed L. Since df(Lt v )/dt = 
vLt"- 1 f'(Lt v ), we must have f'{x) oc x 1 ^- 1 as x ->■ 0, 



leading immediatedly to equation (59). This relation al- 
lows determination of v. 

At K = oo, our model has the same partition func- 
tion as the Heisenberg ferromagnet, and we expect to 
have the same exponents at the transition for K down to 
the multicritical value. Since the magnetization is not a 
gauge-invariant quantity, we instead consider 

9iL > t] = { %V)l L = g{Ln (6o) 

analogous to fi(L,t) (equation 58, with h = 4). We 
use Tr Q 2 here, rather than q 2 (equation 57), because 
the analogous magnetic quantity, m 2 = (./V -1 ^ Si) 2 , is 
nonzero at infinite temperature. Indeed, the presence of 
terms |5i| 2 = 1 in Tr Q 2 , but not in q 2 (recall the discus- 
sion at the beginning of section IX A) result in m 2 ~ 1 /N 
at infinite temperature. Splitting off the i = j terms, 
each of which is \S\ 2 - 1/3 = 2/3, 

(Tr Q 2 ) L = N- 2 • S,) 2 - 1/3) L + 1/N 

= (q 2 )L+N- 1 (l-(q 2 ) L ), (61) 

which tends to 1/N = 1/L 3 as J — ► 0. The nu- 
merical data for this function for L = 8, 12, and 16, 
show intersections indicating a value of the critical cou- 
pling between 0.690 and 0.695, in agreement with the 
best value to date for the three-dimensional Heisenberg 
model [20], J C (K = 0) = 0.693... The shift of the crit- 
ical coupling due to icosahedral anisotropy is evidently 
small. At K = 0.78, the critical coupling is a bit larger: 
J C (K = 0.78) = 0.70. 

For a magnetic system, the finite size scaling of the 
magnetization is also quite useful: 

m = \tff(L/0 = ^"f(Lt"). (62) 

According to eqn. (62), the ratio (5/v is accessible by 
measuring the magnetization, which can be determined 
accurately. Our order parameter is related to the magne- 
tization of an equivalent Heisenberg model, though not 
quite the same. For a finite lattice at values of J not too 
much below J c , the approximation [43] 

((Si ■ S 3 ) 2 ) L = {(S?S?)(S?S?)) L * (S?SP) L (S«SP) L , 

(63) 

should be valid for i and j well-separated. The last ex- 
pression in (63) is proportional to m 4 . Notice that as 
J -s- 0, m 4 ~ 1/7V 2 , but (Tr Q 2 ) ~ 1/N, due to the i = j 
terms in eqn. (61). Thus, (q 2 )h (which lacks these ex- 
tra terms) will behave much more like m 4 than (trQ 2 )L 
does. We conclude that 

(q 2 ) L » LW"f(Lt"). (64) 

At the critical coupling (determined from eqn. 58), 
eqn. (64) can be used to determine ft /v. Combining 
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with v from eqn. (59) gives j3. Repeating this analysis 
at J = 0.693, K = oo, we found v = 0.71 ± 0.04, /3 = 
0.39±0.03, to be compared with the accepted Heisenberg 
values of 0.703 and 0.38 respectively. At K = 0.78, we 
find v = 0.72±0.05 and (3 = 0.4±0.04, the same to within 
error. The scaling in equation (64) can also be checked 
away from the critical coupling, by plotting q 2 L~ if3 l v 
versus L X I V {J — J c ). Tn data collapse is best for the 
Heisenberg values. One expects that effects of a finite 
gauge coupling should be observable at K = 0.78, if they 
are relevant, which they do not appear to be. Scaling 
plots of the order parameter at K = 5 and K = 0.78 are 
shown in figure 9. 

A finite-size scaling analysis was also carried out for the 
specific heat at K = 0.70. The specific heat curves are 
shown in figure 10. Finite-size scaling near a first-order 
transition is quite different from that near a critical point 
[44]. At the discontinuity fixed point governing a first- 
order transition, the only relevant exponent [45] is the 
spatial dimensionality d of the system. One therefore ex- 
pects to see no exponent other than d in the leading-order 
finite size scaling formulas for a first-order transition, in 
contrast to those for a continuous transition which con- 
tain the familiar anomalous exponents. The rounding of 
the delta-function in the specific heat should thus result 
in a peak height linear in the volume of the system, i.e. 
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The same conclusion follows [46] by assuming that config- 
urations of the finite system occur in the partition func- 
tion with a weight which is the average of those appropri- 
ate for the two coexisting phases. The data for K = 0.70 
appearing in figure 10a are consistent with such a scal- 
ing. The data for K = 0.78 (figure 10b), however, do 
not appear consistent with linearity, indicating that the 
transition is not first order at that value of K. 



C. Line Tension 

We also measured the defect line tension in our Monte 
Carlo simulations. The boundary conditions (a) and (b) 
of section V are relatively straightforward to implement 
in the lattice gauge theory of eqn. (11). In two lattice di- 
rections (x and y, say), the boundary conditions are open 
for the spins. The product of the link variables circling 
once around the periphery in the XY plane is required 
to be either +1 (boundary condition a) or —1 (boundary 
condition b). In the remaining lattice directions (just z 
in three dimensions), the boundary conditions are cho- 
sen open with the restriction that either no boundary 
plaquettes are frustrated (a), or only the central plaque- 
ttes on the top and bottom faces are frustrated (b). The 
free energy difference per unit length in the z-direction as 
L x , L y — ► oc is interpreted as a disclination free energy. 
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FIG. 9. Scaling plots for nematic order parameter Tr Q 2 . 
(a) K = 0.78, with v = 0.7, j3 = 0.4. (b) K = 5 (effectively 
infinite), with v = 0.7, /3 = 0.38. The lattice sizes are: L = 8 
(squares), L = 12 (crosses), and L = 16 (circles). 
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Essentially this construction was suggested as a diagnos- 
tic for quark confinement by Mack and Meyer [36] in the 
context of the the Z 2 Higgs-gauge model. 

The general strategy is to force a defect to traverse 
the system in the z-direction, going through the centers 
of the faces at z = and z = L z , but to forbid defects 
from passing through the boundary anywhere else. This 
means fixing the gauge field strength everywhere on the 
boundary. The simplest means of doing this is by actually 
partially fixing the gauge, freezing the link variables on 
the boundary into an appropriate configuration and not 
altering them during the simulation. It is also possible 
to use periodic boundary conditions in the z-direction 
and still force a defect to run through, so that a Wilson 
loop going around the boundary and circling the z-axis 
is —1. The defect is attracted to the boundary by image 
forces, however, where it produces the least disturbance, 
which needs to be avoided. Since the spins do not need 
to be constrained in any way, free boundary conditions 
are used for them. 

The line tension can be studied in such a setup by 
varying L z at fixed L x and L y . As discussed earlier, in 
the infinite volume limit the line tension is expected to 
vanish at K c with the correlation length exponent v\ of 
the Ising model. Finite size effects make interpretation 
of the data difficult. The problem is also exacerbated by 
working on a lattice with boundaries. We measured the 
line tension at J = 0.5 as a function of K. This value 
of J is chosen to be in the nematically disordered region, 
yet with large enough nematic coupling J that its effects 
ought to be observable. The free energy is determined by 
integrating the plaquette expectation value with respect 
to K. In principle we should integrate from the corner 
J = K = 0, but we actually start from K,J= (0,0.5), 
since within statistical error d(SF)/dJ = at that point, 
hence a fortiori for all smaller J. The evolution of the 
line tension with K at J = 0.5 is shown in figure 11. 
There is a kink near the value of K at which the spe- 
cific heat peak occurs. This kink sharpens and moves 
downward with increasing L z , indicating that the ten- 
sion vanishes above K c in the infinite volume limit. The 
derivative of the line tension with respect to J across 
the N/I transition at K = 0.5 is shown in figure 12. As 
the lattice size increases, the peak is sharpening into a 
delta-function reflecting the strictly positive value of the 
nematic order at the transition. Apart from observing 
the qualitative trend with system size, we have done no 
quantitative finite-size scaling analysis. Multiple length 
scales resulting from non-periodic boundary conditions 
makes such an analysis impractical. 



FIG. 10. Evolution of specific heat peaks with lattice size 
in (a) first-order (K = 0.70) and (b) continuous (K = 0.78) 
transition regions. In the first order region, the peak sharpens 
very rapidly with increasing lattice size. The linear scaling of 
peak height with lattice volume is shown in the inset. The 
continuous curves are guides to the eye only. 



D. Defect Line Statistics 

We have also examined the configurations of defect 
lines in our simulations. Figure 13 depicts typical defect 
networks for several representative points in the phase 
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diagram. Upon passing into the totally disordered con- 
finement phase I from the free charge phase T, a dra- 
matic increase in the number of defect loops is evident. 
The increase in total length is not quite so obvious in the 
N/I transition because even in the nematic phase there 
are already many small loops near the transition. When 
the lattice becomes crowded with defective plaquettes, it 
is impossible to unambiguously pick out individual loops 
of defect. The loops coalesce into defect clusters because 
of their intersections. We have compiled some statistics 
for these clusters. Figure 14 shows the distributions of 
total length for the defect clusters in typical configura- 
tions as J is varied at K = 0. The N/I transition is 
characterized not only by a sudden increase in the total 
length of defect, but also by the appearance of very long 
defect loops. The mean square separation of dual lattice 
sites on the clusters was also studied; in a toroidal lattice 
this is the quantity most nearly equivalent to radius of 
gyration. For clusters with radii smaller than the length 
of the lattice, the mean-square separation scales approx- 
imately as net length, as appropriate for a random walk. 
Any self- avoidance of the defect clusters is apparently 
weak. 



FIG. 11. The free energy difference 5F/L between a pair 
of 12 x 12 x L rectangular prisms with and without an ex- 
ternally imposed defect, respectively, as a function of K at 
J = 0.5. Lattice sizes are L=6 (circles), 12 (triangles), and 
16 (squares). Below K ~ 0.7, the line tension is zero in the 
thermodynamic limit. 
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FIG. 12. Derivative of the line tension with respect to J 
across the N/I transition at K = 0.5. The squares are for 
lattice size 10 3 and the triangles for 14 3 . 



E. Higher (Spin)-Dimensional Models 

As indicated at the beginning of this section, we have 
also carried out simulations for the system with four- 
component spins situated at the lattice sites. The phase 
diagram for this case is qualitatively similar, with the or- 
dering transition shifted toward larger J. Note that for 
a four-component nematic point defects are excluded by 
topological considerations (^(HP 3 ) = 0), so that discli- 
nations are the only allowed type of defect. Since there 
is no qualitative distinction, we have not pursued a full 
analysis. 

X. CONCLUSION 

The ordered states of ferromagnetic and nematic media 
are strikingly similar, yet the experimentally observed or- 
dering transitions of the two are quite different. Nematic 
transitions are observed to be weakly first-order and fer- 
romagnetic ones are continuous. 

In this paper, we have shown that the origin of this dif- 
ference lies in the topological differences between the ap- 
propriate order parameter spaces; specifically, the pres- 
ence of disclination lines in the nematic. Furthermore, we 
have shown that this effect can be completely suppressed 
by finite, short-ranged interactions, thereby making the 
disordering transitions of the two systems identical. In 
the process, we have discovered a new phase of nemato- 
genic materials - the phase entered when the nematic 
disorders via a ferromagnet-like, continuous Heisenberg 
transition. This phase exhibits topological order which 
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(b) 




FIG. 13. Snapshot pictures of typical defect configurations, 
(a) in the nematic (JV) phase at (J = 1.3, K = 0.4), (b) in the 
isotropic (i) phase at (J = 1.1, K = 0.4), and (c) in the 
isotropic T phase, at (J = 0A,K = 0.85). Periodicity of the 
boundary conditions is revealed by some of the defects in (a) 
and (b). 
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FIG. 14. Scaling of the radius of gyration (R 2 ) of defect 
clusters with total number of links (N) near the first-order 
transition at K = 0. As J is lowered toward the transition 
value, there is a proliferation of defect loops (both gauge de- 
fects and disclinations) of all lengths. The points fall close to 
the line corresponding to a random walk for all values of J, 
except for some very large clusters, which have wrapped com- 
pletely around the toroidal lattice. The two plots correspond 
to different definitions of defect: (a) is computed from eqn. 
(8) and (b) from the second term of the Hamiltonian in eqn. 
(11). The data (a) was taken with only 20 discrete values of 
the spin because it is necessary that there are no spin values 
with Si ■ Sj = 0; consequently the value of J at the transition 
is somewhat smaller for this case. 
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is destroyed only at a second, distinct transition in the 
Ising universality class. The full scenario is summarized 
in the phase diagram (figure 3). Experimentally observ- 
able critical behavior at the two continuous phase transi- 
tions in that figure are worked out in detail in Paper II. 
The problem of formulating a Landau- Ginzburg descrip- 
tion of our theory and its relation to the usual Landau- 
Ginzburg theory [2] of the nematic will be addressed in 
a future publication. 
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